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We study the stability of singly- and doubly-quantized vortex states of harmonically trapped 
dipolar Bose-Einstein Condensates (BECs) by calculating the low-lying excitations of these conden- 
sates. We map the dynamical stability of these vortices as functions of the dipole-dipole interaction 
strength and trap geometry by finding where their excitations have purely real energy eigenvalues. 
In contrast to BECs with purely contact interactions, we find that dipolar BECs in singly-quantized 
vortex states go unstable to modes with an increasing number of angular and radial nodes for more 
oblate trap aspect ratios, corresponding to local collapse that occurs on a characteristic length scale. 
Additionally, we find that dipolar BECs in doubly-quantized vortex states are unstable to decay 
into a different topological state (with two singly-quantized vortices) for all interaction strengths 
when the trap geometry is sufficiently prolate to make the dipoles attractive, and in windows of 
interaction strength when the trap geometry is sufficiently oblate to make the dipoles repulsive. 



I. INTRODUCTION 

The observation of strong dipolar effects in a Bose- 
Einstein condensate (BEC) of harmonically trapped 
52 Cr 0, Q atoms marks encouraging progress towards 
understanding these novel systems. In contrast to the 
isotropic contact interaction present in condensates of al- 
kali atoms, the dipole-dipole interaction is long ranged, 
anisotropic, and is predicted to induce interesting ground 
state structures [HIH and excitation spectra in both fully 
and partially trapped systems. For example, an exci- 
tation spectrum much like the roton-maxon spectrum 
in superfluid He is predicted in dipolar BECs (DBECs) 
with both three-dimensional (3D) [3j, [5| and quasi-two- 
dimensional (2D) [6( geometries. Additionally, the pres- 
ence of dipolar effects has recently been shown to be crit- 
ical in explaining the rich behaviors of spinor BEC sys- 
tems @, Q j the effects of which are strongly conditioned 
by the attractive part of the dipole-dipole interaction. 
For this reason, it is instructive to compare the dipolar 
system to a BEC with attractive contact interactions. In 
this article, we find the properties of the two systems to 
be in stark contrast. 

To see this contrast, we consider the effects of the 
dipole-dipole interaction on a condensate with a single 
vortex core [j| ■ The conditions for the generation of such 
a DBEC vortex state are studied in Ref. [Tof . First, con- 
sider a trapped BEC with attractive contact interactions, 
characterized by the s— wave scattering length a s of the 
constituent particles and with no vortex. For such a sys- 
tem, there always exists a critical particle number above 
which the condensate goes unstable, with preference to 
collapse in the region of maximum density at the cen- 
ter of the trap [ll|, [l2| ■ Stirring the condensate into a 
vortex state serves to stabilize the system by introducing 
a kinetic energy component due to angular momemtum 
that offsets the interparticle attraction. So, in general, 
the vortex will sustain a larger number of particles than 
the non-vortex state, and is more stable. 

The case for a DBEC, however, is quite different. 
Consider a DBEC with its dipolar entities polarized in 



the z-direction. Because the dipole-dipole interaction 
is anisotropic, the structure and stability of a DBEC 
depends strongly on the trap aspect ratio A = w 2 /w p , 
where lo z and u p are the axial and radial trap frequencies, 
respectively. For smaller A, the dipole-dipole interac- 
tion distends the condensate into a prolate shape, where 
macroscopic collapse can occur due to long-range attrac- 
tion in the direction of polarization. As A is increased, 
the condensate is stabilized since the dipole-dipole inter- 
action is predominitely repulsive in this more oblate ge- 
ometry. For a moderate number of particles, the energy 
cost of stacking the dipoles in the direction of polariza- 
tion is outweighed by the tight harmonic confinement in 
this direction. However, for a sufficiently large number 
of particles, a DBEC in an oblate trap is subject to an 
instability due to local density fluctuations, which are for- 
eign to the contact potential case. The attractive part of 
the dipole-dipole interaction dominates in regions where 
the higher-density fluctuations occur, inititating a local 
collapse of the condensate. As will be discussed below, 
this instability is intimitely connected with an excitation 
that goes soft at a critical number of dipoles, and that 
has been dubbed the "discrete-roton." Q Signatures of 
this local collapse have been articulated in the simula- 
tions of Ref. [13|] , where the manifestation of the roton in 
the collapse of a DBEC is shown for the early stages of 
collapse. 

In the presence of a singly-quantized vortex, the region 
of high density is forced away from the center of the trap 
due to the zero-density of the vortex core. Depending 
on the aspect ratio A of the trap, this either serves to 
stabilize (for smaller A) or destabilize (for larger A) the 
DBEC. For smaller A, the vortex core simply breaks the 
prolate shape of the condensate along the direction of 
polarization, eliminating much of the attractive dipole- 
dipole interaction in this direction and thus increasing 
the energy due to interactions. Conversely, for larger A 
the vortex increases the density in the periphery of the 
core and thus encourages local collapse. Just as the ro- 
ton wavelength is set by the the confinement length in 
the direction of polarization (z), the local density flue- 
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tuations occur at the same length scale regardless of the 
trap geometry. Widening the trap radially while keeping 
the axial trapping frequency fixed makes more room for 
regions of density fluctuations instead of enlarging the ex- 
isting regions. This marks a clear and important distinc- 
tion between the dipole-dipole and contact interactions. 
An additional consideration relevant to the stability of 
DBECs with a vortex are the excitations of the vortex 
core itself [14|. As we will see below, these excitations 
are unlikely to play a role in destabilizing the vortex for 
oblate, or even mildly prolate, traps. 

Our paper is organized as follows: In section [II] we de- 
scribe the model that we use to study the ground state 
with a vortex, including the novel algorithm that we em- 
ploy to carry out our calculations. In section IIIII we 
discuss the Bogoliubov de Gennes (BdG) formalism and 
present our calculations of the BdG spectrum in reference 
to the stability of the system. Finally, in section HVl we 
calculate the BdG spectrum of a DBEC with a doubly- 
quantized vortex. Interestingly, in a more prolate trap 
geometry where the dipoles are mostly attractive we find 
that the BdG spectrum looks very similar to the BEC 
with purely attractive contact interactions, while in a 
more oblate trap geometry where the dipoles are mostly 
repulsive we find that the BdG spectrum looks very sim- 
ilar to that of the BEC with purely repulsive contact 
interactions, having windows of dynamical stability for 
certain dipole-dipole interaction strengths [HI, [l6[ . 



II. METHODS 

At ultracold temperatures, a condensate of N bosons 
trapped in an external potential U(r) may be described 
within mean field theory [l7j by the nonlocal Gross- 
Pitaevskii Equation (GPE): 



x J dr'%(r')V(r-r')q {r') 



*o(r) = M*o(r) (1) 



where Vfol 1 ") is the condensate wave function with unit 
norm, r is the distance from the trap center, M is the bo- 
son's mass and V(r— r') is the two-particle interaction po- 
tential. For a cylindrically symmetric harmonic trap, the 
external potential is given by U(r) — \Mw 2 {p 2 + X 2 z 2 ), 
where A = uj z /uj p is the trap aspect ratio. The interaction 
potential for a dipolar system has the form [l8[ 



V(r-r') = 
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where d is the dipole moment and is the angle between 
the vector r — r' and the dipole axis. To illuminate purely 
dipolar effects, we set a s = in this work, a limit that has 
been achieved experimentally in a BEC of atomic 52 Cr 
by employing a Fano-Feshbach resonance |l9j. 



To characterize the structure and stability of a DBEC 
with a vortex, we introduce the dimensionless character- 
istic dipole strength, 



n a ho 



(3) 



where aha — \/Ti/Mlu p is the radial harmonic oscilla- 
tor length. Notice that increasing (decreasing) D corre- 
sponds to either increasing (decreasing) the number of 
particles in the condensate or the square of the dipole 
moment of the particles. So for a DBEC of, say, 52 Cr, 
one must change the number of atoms in the condensate 
in order to change D, since the magnetic dipole moment 
of 52 Cr (6 hb) is effectively fixed. 

The second term in Eq. ([2]) describes the two-body 
dipole-dipole interaction for dipoles that are polarized 
along the trap axis (z-axis) [20j, as may be achieved 
by applying a strong external field to the condensate. 
This term is long ranged (oc 1/V 3 ), anisotropic, and 
has a coordinate-space singularity as the distance be- 
tween the particles goes to zero. These concerns are 
handled by treating the mean field dipole term in the 
GPE, d 2 J dr'frg(r') 1 -^£^ V (r'), by means of convo- 
lution. This method eliminates the singularity in coor- 
dinate space since the dipolar mean field in momentum 
space, Vb(k), is given by [2ll ] 



~ 4-7T 

V D (k) = Y (3kl/k 2 -l) 



(4) 



The dipolar mean field in coordinate space may then be 
calculated in terms of Vb(k), 



<f> D (r)=T- 1 Vb(k)n(k) 



(5) 



where J- 1 is the inverse Fourier transform operator and 
n(k) is the Fourier transform of the condensate density. 

In general, these transforms must be computed in three 
dimensions to capture the three-dimensional (3D) nature 
of the system. However, the system that we are consider- 
ing possesses cylindrical symmetry in both the trapping 
and dipolar interaction potentials. With such a symme- 
try, the condensate states may be written in cylindrical 
coordinates as eigenstates of the orbital angular momen- 
tum projection s, *o( r ) = i>o{p, z)e lsv [22|]. The s = 
state corresponds to a rotationless condensate, while the 
s > states correspond to condensates with vortices of 
charge s, or equivalently, with hs units of orbital angu- 
lar momentum per particle. This formulation allows for 
a distinctive computational algorithm to be applied to 
the problem, reducing a fully 3D calculation to a 2D one 
by working in cylindrical coordinates and integrating out 
the simple ip dependence of the state. Specifically, the 
algorithm uses a one-dimensional (ID) Hankel transform 
of order s in the p-coordinate and a ID Fourier trans- 
form in the z-coordinate to transform a function with 
the angular dependence e lsv into momentum space. For 
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example, the transform of the wave function ipaiP, z)e lsip 
is 



dze %kzZ 



V>o(k) = V2irr s e 
pdpi/j (p,z)J s (k p p) 



(6) 



where J s {k p p) is the Bessel function of order s. In prac- 
tice, the Fourier transform and the Hankel transform of 
order s are performed on discrete grids. The applica- 
tion of this algorithm to solving the GPE is detailed in 
Ref. [13 . 

In addition to calculating the ground state with a spec- 
ified vorticity using this algorithm, we also extended it to 
calculate the BdG excitation spectrum in the presence of 
a vortex ground state (see section IIIII) , which describes 
the low-lying excitations and reveals critical information 
regarding the dynamical stability of the system. This 
extension of our algorithm is described in Appendix [AJ 



III. EXCITED STATES AND STABILITY 



To study the elementary excitation spectrum of a 
DBEC with a given projection of angular momentum s 
in the ground state, we use the BdG ansatz and write 
a wave function of the time-dependent GPE of the form 
(with h = 1) 



(7) 



where p is the chemical potential of the ground state, 
ipo(p, z)e lsip , and $(r,£) is given by 

(r, t) = S \u(p, z y( m v^t) + v * ^ z ) e -<(mv.-«t)l ( 8 ) 

where 6 <C 1 to ensure that the excitations have small 
amplitudes, uj is the frequency of the excitation and 
the modes u(p, z)e lmv and v*(p, z)e~ lmv are normalized 
by |H 



dr[\u(p,z)\ 2 -\v(p,z)\ 2 ] =1. 



(9) 



The BdG modes are characterized by the quantum num- 
ber m, being their projection of orbital angular momen- 
tum onto the z-axis. This ansatz represents a vortex 
with angular momentum Tis per particle giving rise to 
excitations with angular momentum h(s ± to). By in- 
serting Eq. J7]) into the time-dependent GPE (Eq. ([T]) 
with p on the right hand side replaced by and lin- 

earizing about S, the coupled BdG equations are derived 
by collecting terms evolving in time like e~ iU}t and e ZUJt , 
respectively, 



Lju{p,z)e lmv = [Hq — p 
drV*(p', z')V N (r - v')Mp', A z> imv 

dr>r (p', z')V N (v ~ v')u(p', z')e m ^'MP, *) 

dr'v(p>, z')V N {v - t')Mp\ z')e m *'MP, z) 

-Ljv(p,z)e im * = [H -p 



(10) 



+ / dv'r {p',z')V N {v-v')Mp'^') 



v(p,z)e l 



dr'Mp', z')V N (v - r')v(p', z')e lm * > *(p, z) 
dr'u(p\ z')V N (r - r')r (p\ z')e^ >*(p, z), 



(11) 



where H = -^gV 2 + U(r) and V N (r - r') = (N - 
l)V r (r — r'). By using Hankel transforms to compute 
the interaction terms in Eqs. (| 10[) and (fTTj) . we are able 
to account for the angular dependence of the integrands 
by using a Bessel function of the appropriate order, as 
described in Appendix [5] To solve these equations, we 
write them in matrix form, as in fl5l . [23| , and diagonalizc 
them numerically to find the eigenvectors (ui, v*) and the 
eigenvalues Wj. 

We point out that while there are solutions of the BdG 
equations of the form (tij, u*), there are always solutions 
of the form (v*,Ui) with Ui —¥ —u>i and m —¥ —m. 
For the case of s = BECs, there is a solution of 
the original form (ui,v*) with the same uji but with 
m — > —m. This simply expresses the fact that counter- 
rotating excitations are degenerate due to the reflection 
symmetry of the s = ground state. The presence of 
a vortex breaks this double degeneracy. We shall say 
that the excitation (ui,v*) has a positive norm when 
J dr [\u(p, z)\ 2 — \v(p, z)| 2 ] > 0. It can then be normal- 
ized such that J dr [\u(p, z)\ 2 — \v(p, z)| 2 ] =1. The solu- 
tion (v* , Ui) with u>i — > —uii and m — > —to will then have 
a negative norm, obeying J dr [\u(p, z)\ 2 — \v(p, z)| 2 ] = 
— 1. A positive norm mode with a negative energy eigen- 
value signifies that there exists a lower energy solution of 
the GPE; the same situation is represented by a negative 
norm mode with a positive energy eigenvalue [24j . 

The solutions of the BdG equations characterize the 
stability of s = 1 DBECs. The global thermodynamical 
instability of s = 1 DBECs is seen as a negative norm 
BdG mode with to = 1 and positive energy for all trap as- 
pect ratios and dipolar interaction strengths. This mode 
corresponds to the system's decay into the energetically 
favored rotationless ground state, just as for BECs with 
purely contact interactions. The component of the mode 
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FIG. 1: (color online). The black dashed line marks the maxi- 
mum dipole strength D, for a given trap aspect ratio A, above 
which the Gross-Pitaevskii Energy functional has no mini- 
mum corresponding to the s = 1 DBEC. The solid line marks 
a more restrictive stability line, determined by the onset of 
dynamical instability, signaled by the emergence of an imagi- 
nary energy in an excitation mode. The pink (darker) region 
represents where oscillations with local minima are observed 
on dynamically stable states. The inset is an isodensity sur- 
face plot of an s = 1 DBEC at the point in parameter space 
indicated by the arrow. The ripples in the density are ex- 
plained in Ref pj. 



with angular dependence e - l ( TO - s )y — l is in this case ro- 
tationless, capturing the symmetry of the s — ground 
state. Since this mode populates the core of the vortex, 
it is rcfcred to as a core mode. However, at ultracold 
temperatures, thermodynamical stability is less relevant 
in characterizing the stability of a condensate since there 
needs to be some thermal processes acting on the system 
to dissipitavely drive it into a lower energy state. We 
therefore disregard thermodynamical instability in the 
following. 

Instead, we focus on studying the dynamical stabil- 
ity of s = 1 DBECs. The emergence of a BdG energy 
eigenvalue with a nonzero imaginary part corresponds 
to a dynamical instability in the system [22| . For ex- 
ample, suppose that a BdG mode (iii,v*) has energy 
LUi = ujr + Wi with u>i 0; then the mode will have 
the time dependences e ("i—«"ii)t and e - either 
exponentially growing or decaying in time. Consequently, 
lui determines the rate of decay of the condensate, given 
by r = 

We determine where s = 1 DBECs are dynamically 
stable by finding the region in parameter space where all 
of the BdG modes have purely real energy eigenvalues. 
This region is shown by the colored portion of Figure Q] 
The dashed line in this figure marks, for a given A, the D 
below which we find a local minimum of the CP energy 
functional [23| by using our reduced 2D algorithm. We 
find that, for all A, the D above which the GP energy 



functional has no minimum corresponding to an s = 1 
ground state and the D at which the BdG spectrum be- 
gins to possess imaginary energy, denoted D C rit, are never 
equal. Indeed, dynamical instability occurs for values of 
D at which the GPE has a solution. This is because 
in the 2D minimization of the vortex-state energy, per- 
turbations that break the s = 1 angular momentum are 
not allowed, these are only examined later with the BdG 
equations. Using a fully 3D calculation, we check the 
accuracy of D C rit for various trap aspect ratios by time 
evolving the condensate wave function with an initial ran- 
dom perturbation. The -D C rit that we calculate using the 
3D algorithm, corresponding to the D at which we ob- 
serve collapse under time evolution, agrees with the D CI n 
that we calculate by finding imaginary energy eigenvalues 
in the BdG spectrum using our 2D algorithm. The pink 
(darker) region in Figure [T] represents the region where 
we find dynamically stable s — 1 ground states having ra- 
dial ripples with local minima, as illustrated by the inset. 
This feature has been explained in detail in Ref. @. 

We find that s — 1 DBECs possess imaginary energy 
in their BdG spectrum only when two modes of oppo- 
site norm are degenerate with each other, just as is the 
case for a BEC with contact interactions. This circum- 
stance was recently studied in Ref. [25| . where it is con- 
firmed perturbatively for BECs with contact interactions. 
Ref. 15J also confirmed this claim using a two-mode ap- 
proximation. Indeed, we find that the same holds true 
for DBECs, where the only difference between the two 
systems is the shape of the mean field potential. At all 
aspect ratios, we observe, for some finite value of D, two 
modes with opposite norm approach and then go degen- 
erate with each other at _D cr it . At the point of degeneracy, 
the modes develop equal and opposite imaginary ener- 
gies, signifying dynamical instability. If two modes that 
have the same norm approach each other, they undergo 
an avoided crossing instead of becoming degenerate. 

For a BEC with pure contact interactions in the s = 1 
vortex state, the mode that defines the onset of dynami- 
cal instability is independent of aspect ratio A. Positive 
contact interactions ensure dynamical stability while neg- 
ative contact interactions (for A > 0.3) bring about a 
dynamical instability due to an m = 2 mode [26|. Addi- 
tionally, the s = 2 state is dynamically unstable due to an 
ttl — 2 mode for negative contact interactions, while an 
m = 2 mode defines windows of dynamical stability for 
positive contact interactions. This holds true for these 
systems no matter how oblate the trap. 

The case for a DBEC, however, is quite different. Fig- 
ure[5]illustrates the imaginary parts of the BdG spectrum 
for m = 1 — 5 for DBECs in traps with aspect ratios 
A = 2 and A = 15. Where these imaginary energies are 
zero, from D = to -D C rit, the condensates are dynam- 
ically stable. Notice that for A = 2, an m = 2 mode 
develops imaginary energy at a I? well below the other 
modes, defining £> C rit for this aspect ratio. However, an 
rn = 4 mode serves to define £> cl -it for A = 15. 

Indeed, unlike BECs with contact interactions, modes 
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FIG. 2: (color online), a) The imaginary part of the BdG 
excitation spectrum for a number of m values for a DBEC 
in a trap with aspect ratio a) A = 2 and b) A = 15. For 
A = 2, the m — 2 modes clearly develop imaginary energy 
eigenvalues at a D smaller than the other modes, defining 
D cr it ~ 8 for this aspect ratio. For A = 15, the m = 4 
modes develop imaginary energy at a D smaller than the other 
modes, defining D CI n ~ 108 for this aspect ratio. 



with different m quantum numbers serve to define -D cr it 
at different aspect ratios for DBECs. For moderate trap 
aspect ratios (such as A = 2), an m = 2 mode defines 
£> cr it for the DBEC, similar to the case for contact inter- 
actions. However, as the trap aspect ratio is increased 
to more oblate shapes, modes with larger m quantum 
numbers develop imaginary energy eigenvalues at smaller 
values of D than the m = 2 mode. Figure [21a) illus- 
trates this by plotting the differences between the £)'s at 
which the BdG modes with different angular symmetries 
first develop imaginary energy eigenvalues, and D CIlt , as a 
function of A. Thus, for a given A the lowest curve identi- 
fies the symmetry of the unstable mode. For 6 < A < 12, 
anm = 3 mode defines D CT i t while for larger aspect ra- 
tios, an m — 4 mode defines -D C rit- Although it is not 
shown here, we find that at even larger aspect ratios the 
vortex decays into still higher m-modes. 

The relevance of the m-dependent dynamical instabil- 
ity is that the dipole-dipole interaction leads a BEC to 
instability locally and at a fixed length scale, the wave- 
length of which is determined by the axial harmonic os- 
cillator length. We find that, at the onset of imaginary 
energy, these modes have radial nodal spacings very simi- 
lar to that of the roton o n the rot ationless DBEC, namely 
A/2 ~ 7rZ z , where l z = \J h /M ui z is the axial harmonic os- 
cillator length. The angular dependence of these modes 
behaves in the same way. Increasing A decreases the ra- 
tio l z /ah , so more radial nodes, fixed by l z , can fit into 
the condensate for larger A. In the same way, more an- 
gular nodes can fit into the condensate, therefore bring- 
ing about dynamical instability by modes with larger m 
quantum number, and hence more angular nodes. 



FIG. 3: (color online), a) The difference between the D at 
which BdG modes with different m quantum numbers first de- 
velop imaginary energy eigenvalues and D cr it, the smallest D 
at which any mode develops an imaginary energy eigenvalue. 
Modes with larger m serve to define D cr it for more oblate 
traps, b) The same, but for negative contact interactions in- 
stead of dipole-dipole interactions. An m — 2 BdG mode is 
always the first to develop an imaginary energy eigenvalue for 
this case, except for trap aspect ratios A < 0.3, for which an 
m = 1 mode plays this role. 



All of the previously discussed BdG modes that we 
identify as being responsible for dynamical instability 
are axially symmetric and nodeless in z. Modes that 
break this axial symmetry can correspond to vortex ex- 
citations, where the vortex core itself may tilt or bend, 
and have been termed "kelvon" modes. Ref. [l4| reports 
that, for a singly-quantized vortex in a DBEC that is 
otherwise spatially homogeneous, the condensate is dy- 
namically unstable to a kelvon mode when an external 
periodic potential is applied along the direction of the 
vortex. We find that, in a harmonically trapped DBEC, 
a mode with a single node at z = determines D c rit for 
A < 0.28. Modes of this type might therefore correspond 
to a kelvon-instability in prolate traps, but we leave these 
considerations for future work. 

As was done in Ref. [26[ for self-attractive BECs 
in the singly-quantized vortex state, we perform time- 
dependent simulations of a DBEC where D is chosen to 
be just above -D cr it, enabling us to go beyond the small 
deviations from the stationary vortex state and see the 
actual process of collapse. Initializing the simulations 
with random noise, we observe collapse, at all aspect ra- 
tios, with an angular symmetry corresponding to the m 
quantum number of the mode that first develops an imag- 
inary energy eigenvalue. 
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IV. STABILITY OF VORTICES WITH s = 2 

The dynamical instability of condensates with doubly- 
quantized vortices and contact interactions has been 
studied extensively [TH HH, H3, HH . These studies report 
windows of positive scattering length where the BECs 
are dynamically unstable to an m — 2 BdG core mode, 
as well as dynamical instability for all values of negative 
scattering length due to an m = 2 mode. Knowing that 
the dipolar mean held in a DBEC can be engineered to 
be more attractive or repulsive for smaller or larger trap 
aspect ratios, respectively, we investigate the presence of 
these features in s = 2 DBECs. When the harmonic trap 
is more spherical, the dipoles are free to stack vertically, 
creating an attractive mean field in this direction. How- 
ever, in pancake shaped traps the dipoles create a more 
repulsive mean field. Thus, for larger trap aspect ratios 
DBECs are more self-repulsive than for smaller aspect 
ratios, mimicking the mean field of condensates with re- 
pulsive contact interactions. As an example, we calculate 
the contribution of the dipolar mean field to the energy 
of a DBEC in a trap with aspect ratio A = 2 and with 
A = 15 for a fixed D = 10. In the A = 15 trap, we find 
that this contribution is about five times larger than in 
the A = 2 trap. 

Indeed, for a DBEC in a trap with aspect ratio A = 2 
we find that there exists an m = 2 mode with a com- 
plex energy eigenvalue for all values of D. However, for 
A = 15 we find that there are windows in D where an 
to = 2 BdG mode has a complex energy eigenvalue, while 
this same mode has purely real energy outside of these 
windows, as illustrated in Figure |4] For trap aspect ra- 
tios A < 7.5, there are no windows of dynamical stabil- 
ity and the condensate is dynamically unstable for all 
D. However, windows of dynamical stability appear for 
aspect ratios A ^ 7.5 and continue for larger A. As is 
reported in Ref. [15[ , we find that there is an to = 2 core 
mode with negative norm and positive real energy that 
increases monotonically as it goes successively degener- 
ate with positive norm modes as D is increased to create 
the windows of dynamical instability. This mode repre- 
sents the s = 2 condensate's instability to splitting into a 
condensate with two singly quantized vortices. The core 
mode is thermodynamically unstable for all values of D 
and is only dynamically unstable for the windows shown 
in Figure |4j 
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FIG. 4: (color online). The positive imaginary part of the 
energy eigenvalues for a doubly quantized (s = 2) DBEC in 
a trap with aspect ratio a) A = 2 and b) A = 15. For A = 2, 
an m = 2 mode possesses imaginary energy for all values of 
D, signifying a dynamical instability for all D. However, for 
A = 15 there are windows of dynamical stability as the m — 2 
mode alternates having and not having imaginary energy for 
different values of D, as is the case for a purely repulsive s = 2 
BEC. 



complex energy eigenvalue. The value of D at which this 
imaginary energy appears marks the threshold of dynam- 
ical stability, -D C rit, for the given trap aspect ratio. By 
inspecting the BdG spectrum for all to quantum num- 
bers, we determine the symmetry of the mode that is 
responsible for the dynamical instability in the conden- 
sate. We find, in contrast to BECs with purely contact 
interactions, that DBECs with a singly-quantized vortex 
go unstable to modes with larger to quantum numbers 
for larger trap aspect ratios, signifying a type of local col- 
lapse of these condensates. We have checked the accuracy 
of D CT n for various trap aspect ratios by performing fully 
3D simulations. Additionally, we find that, in analogy 
to a self-repulsive BEC with a doubly-quantized vortex, 
at larger trap aspect ratios there are successive regions 
in D where the s = 2 DBECs are dynamically unstable 
due to an to = 2 core mode, while the condensates are 
dynamically stable outside of these regions. 
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Appendix A: Hankel Transforms and Interpolation 

Consider the Hankel transform of a function f(p), f(k), 
where f(p) ~ for p > R and f(k) ~ for k > K, and 
define S = RK. The discrete Hankel transform (DHT) 
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of order to of f(p) is then given by |2S 

N 



f(k m i) 



2_ \ " fiPmj) 



-Jr, 



5' 



(Al) 



where i,j = 1,...,N, p mj = a mj /K, k ml = a ml /R and 
a m i is the i th root of J m (p). Conversely, the inverse DHT 
of order to of f{k) is given by 



2 " 



m.) > 



? = 1 Jrn+l( a mj) 



( amj J" M ) (A2) 



Eqs. (|A1[) and (|A2|) show that in order to perform a DHT 
of order m on a function, the function should be defined 
on a grid proportional to the zeros of the m-order Bessel 
function, p m i. 

For the problem we are considering, the functions 
f(p) and f(k) have the angular dependence e lmv and 
e -imk, p ^ respectively, which is why the DHTs above are 
written with Bessel functions of order m. When we 
calculate the BdG modes for the vortex states of a 
DBEC (see section ITTT|) . we take DHTs of functions like 
ipo(p, z)e~ lsv u(p, z )e' l( - m+s ' >v where i/'oG ' z) is defined on 
a grid proportional to the zeros of the s-order Bessel func- 
tion and u{p, z) is defined on a grid proportional the the 
zeros of the Bessel function of order to + s. To perform a 
DHT on the p-coordinate of the function ifio(p, z)u(p, z), 
the function must be defined on a grid proportional to 
the zeros of the Bessel function of order to, so ^Jq(p,z) 
and u(p, z) must both be interpolated onto this grid. To 
accomplish this, we have developed an accurate interpo- 
lation scheme based on the DHT itself. 

As explained in Ref. [3(j, the function f(p) may be 
expanded in an m th order Bessel series, 



f(p) 



N 

£< 

i=l 



i J n 



where the coefficients c m i are given by 



: 7 via / f{p)Jm (a ml ^) pdp. 

m+l{ami)\ Jo V 



" ll i? 2 [J, 

(A4) 

The integral in Eq. (|A4j) is just the Hankel transform ([6]) 
with a m i/R — k m i, giving the transformed function 
f(ki). If p is discretized in Eq. (|A3|) . then this prescrip- 
tion gives exactly Eq. (|A2jl . 



We wish to consider the case where our function is 
defined on the grid p™, proportional to the zeros of the 
Bessel function of order m, but it needs to be defined on 
the grid p n i, with n =/= to. To do this, we expand f(p n i) 
in a Bessel series, 



f(Pni) 



N 



( Pru_ 
[Omi R 



(A5) 



where the coefficients c m j are given by Eq. (|A4[) and are 



computed in terms of the zeros of the Bessel function 
of order to. However, in Eq. (|A5|) . the function is ex- 
panded on the grid p ni , proportional to the zeros of the 
Bessel function of order n. The interpolation algorithm 
then simply follows by inserting the expression for the 
coefficients, 



2 ^ 

f(Pni) = ~fy2 



mj ) 



7 - =1 Jm+l( a ™j) 



where f(k m j) is the discrete Hankel transform of f(p m j). 



(A3) Eq. (TAll . 
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